BeginPackage["CoherentCluster`",{"MathMR`","MRData`","SpinBath`","NVInteraction`"}];

Bx					::usage="Effective field in x direction of nuclear spin.";
By					::usage="Effective field in y direction of nuclear spin.";
Bz					::usage="Effective field in z direction of nuclear spin.";
xx					::usage="Dipolar interaction component xx";
yy					::usage="Dipolar interaction component yy";
zz					::usage="Dipolar interaction component zz";
xy					::usage="Dipolar interaction component xy";
yz					::usage="Dipolar interaction component yz";
zx					::usage="Dipolar interaction component zx";
yx					::usage="Dipolar interaction component yx";
zy					::usage="Dipolar interaction component zy";
xz					::usage="Dipolar interaction component xz";

Spin2Cluster		::usage="Calculate the coordinates of a 2-spin cluster of given index.";
Spin3Cluster		::usgae="Calculate the coordinates of a 3-spin cluster of given index.";
CoordConvert		::usage="";
ClusterInfo			::usage="Show the cluster information.";
ClusterInfoCoord	::usage="Show the cluster information.";

MediatedTensor		::usage="";
HamiNSpinXYZ		::usage="Generate the Hamiltonian form for given spins.";
NumericClusterHami	::usage="Calculate the numerical value of a Hamiltonian for given hamiform and interaction strength.";
ClusterHami			::usage="Calculate the numerical value of a Hamiltonian for given bath spin coordinates and central spin state.";
NoiseMat			::usage="Calculate the numerical value of a noise matrix for given bath spin coordinates and central spin state.";
BathMat				::usage="";

TimeIntervalSeq		::usage="Calculate the contour time intervals for the given pulse timings.";
FilterFunc			::usage="Calculate the filter function of given pulse seq..";
FilterFunction		::usage="";
FilterFunctionArray	::usage="";

PostFunction		::usage="";
ClusterCoherence	::usage="Calculate the cluster coherence for a given pulse seq., conditional Hamiltonians, and evolution time list.";
ClusterEvolution	::usage="ML realization of cluster evolution";
CoherentClusterEvolution			::usage="Calculate the cluster coherence for given cluster (index or coordinates), pulse seq., and evolution time list.";
CoherentClusterEvolutionGaussian1	::usage="";

NoiseAnalysis		::usage="";
NoiseFromState		::usage="";
CoherentClusterEvolutionGaussian	::usage="";

ReadCCEData			::usage="Read .cce file generated by CCE programe.";
PlotCCECluster		::usage="Plot cluster evolutions stored in .cce file.";
GlobalCCEInfo		::usage="Show the informations of .cce data.";

(*=================================================================*)
(*Cluster Coordinates*)
Begin["`ClusterCoord`"];

CoordConvert[coord1_]:=
	Module[{res,element,spinvalList,gammalist,coord},
		element=coord1[[All,1]];
		spinvalList= (element/.ElementRule)[[All,1]];
		gammalist= (element/.ElementRule)[[All,2]];
		coord=coord1[[All,{2,3,4}]];
		res={coord,spinvalList,gammalist};
		Return[res];
	]
	
Spin2Cluster[{i_,j_,k_},r_,d_]:=
	Module[{res,baseCoord,shift,res1,res2},
		baseCoord={i,j,k};
		shift=$DiamondPairList1[[r,d]];
		res1=baseCoord+shift[[1]];
		res2=baseCoord+shift[[2]];
		res={res1,res2}*$DiamondLatticeConst;
		Return[res];
	]

Spin3Cluster[{i_,j_,k_},r_,d_]:=
	Module[{res,baseCoord,shift,res1,res2,res3},
		baseCoord={i,j,k};
		shift=$DiamondTriSpinList[[r,d]];
		res1=baseCoord+shift[[1]];
		res2=baseCoord+shift[[2]];
		res3=baseCoord+shift[[3]];
		
		res={res1,res2,res3}*$DiamondLatticeConst;
		Return[res];
	]

ClusterInfo[nspin_,{{i_,j_,k_},r_,d_},state_,{b_,theta_,phi_},gammaBath_:$GammaC13]:=
	Module[{res,res1,res2,coord,pairs,hf,dip,dist,orient,hfdiff,gammaRatio},
		gammaRatio=gammaBath/$GammaC13;
		
		coord=If[nspin==2,Spin2Cluster[{i,j,k},r,d],Spin3Cluster[{i,j,k},r,d]];
		pairs=Subsets[coord,{2}];
				
		dist=Norm[Mean[coord]-$QubitPos];
		orient={"Orientation"}~Join~(Flatten[Differences[#]]&/@pairs);
		hf=gammaRatio*HyperfineField[#,state,b,theta,phi,3]&/@coord;
		
		dip={"dip"}~Join~((gammaRatio^2*DipolarInteraction[#[[1]],#[[2]],b,theta,phi,3])&/@pairs);
		hfdiff={"HF diff"}~Join~(Flatten[Differences/@Map[(gammaRatio*HyperfineField[#,state,b,theta,phi,3])&,pairs,{2}]]);
		
		res1=Grid[{{"dist",dist}~Join~Table[SpanFromLeft,{nspin-1}],{"Coord"}~Join~coord,{"HF"}~Join~hf},Frame->All];
		res2=Grid[{{"index"}~Join~Subsets[Range[nspin],{2}],orient,hfdiff,dip},Frame->All];
		
		res={res1,res2};
		Return[res];
	]/;(nspin==2||nspin==3)

ClusterInfo[nspin_,coord_,state_,{b_,theta_,phi_},gammaBath_:$GammaC13]:=
	Module[{res,res1,res2,pairs,hf,dip,dist,orient,hfdiff,gammaRatio},
		gammaRatio=gammaBath/$GammaC13;
		
		pairs=Subsets[coord,{2}];
		
		dist=Norm[Mean[coord]-$QubitPos];
		orient={"Orientation"}~Join~(Flatten[Differences[#]]&/@pairs);
		hf=(gammaRatio*HyperfineField[#,state,b,theta,phi,3])&/@coord;

		dip={"dip"}~Join~((gammaRatio^2*DipolarInteraction[#[[1]],#[[2]],b,theta,phi,3])&/@pairs);
		hfdiff={"HF diff"}~Join~(Flatten[Differences/@Map[(gammaRatio*HyperfineField[#,state,b,theta,phi,3])&,pairs,{2}]]);

		res1=Grid[{{"dist",dist}~Join~Table[SpanFromLeft,{nspin-1}],{"Coord"}~Join~coord,{"HF"}~Join~hf},Frame->All];
		res2=Grid[{{"index"}~Join~Subsets[Range[nspin],{2}],orient,hfdiff,dip},Frame->All];

		res={res1,res2};
		Return[res];
	]/;nspin>=4

ClusterInfoCoord[nspin_,coord_,state_,{b_,theta_,phi_},gammaBath_:$GammaC13]:=
	Module[{res,res1,res2,pairs,hf,dip,dist,orient,hfdiff,gammaRatio},
		gammaRatio=gammaBath/$GammaC13;
		
		pairs=Subsets[coord,{2}];
		
		dist=Norm[Mean[coord]-$QubitPos];
		orient={"Orientation"}~Join~(Flatten[Differences[#]]&/@pairs);
		hf=(gammaRatio*HyperfineField[#,state,b,theta,phi,3])&/@coord;

		dip={"dip"}~Join~((gammaRatio^2*DipolarInteraction[#[[1]],#[[2]],b,theta,phi,3])&/@pairs);
		hfdiff={"HF diff"}~Join~(Flatten[Differences/@Map[(gammaRatio*HyperfineField[#,state,b,theta,phi,3])&,pairs,{2}]]);

		res1=Grid[{{"dist",dist}~Join~Table[SpanFromLeft,{nspin-1}],{"Coord"}~Join~coord,{"HF"}~Join~hf},Frame->All];
		res2=Grid[{{"index"}~Join~Subsets[Range[nspin],{2}],orient,hfdiff,dip},Frame->All];

		res={res1,res2};
		Return[res];
	]

End[];

(*=================================================================*)
(*Cluster Hamiltonian*)
Begin["`ClusterHami`"];

HamiNSpinXYZ[ncl_,Jvallist_]:=
	Module[{res,fullId,Idxlist,pairIdxlist,lpair,
			singlebodyZ,singlebodyX,singlebodyY,singleHamiZ,singleHamiX,singleHamiY,singleHami,
			xxRules,xxCoupling,xxHami,yyRules,yyCoupling,yyHami,zzRules,zzCoupling,zzHami,
			xyRules,xyCoupling,xyHami,yxRules,yxCoupling,yxHami,yzRules,yzCoupling,yzHami,
			zyRules,zyCoupling,zyHami,zxRules,zxCoupling,zxHami,xzRules,xzCoupling,xzHami,
			xOpList, yOpList, zOpList, idOpList, dimList},
		dimList=(2#+1)&/@Jvallist;
		xOpList=SpinOperatorX/@Jvallist;
		yOpList=SpinOperatorY/@Jvallist;
		zOpList=SpinOperatorZ/@Jvallist;
		idOpList=IdentityMatrix/@dimList;
		
		If[ncl==1,
			res=Subscript[Bz, 1]zOpList[[1]]+Subscript[Bx, 1]xOpList[[1]]+Subscript[By, 1]yOpList[[1]];
			,
			(*ncl > 1 *)
			Idxlist=Range[ncl];
			pairIdxlist=Subsets[Idxlist,{2}];
			lpair=Length[pairIdxlist];
			fullId=idOpList;(*Table[Id,{ncl}];*)
			
			
			(*single body hamiltonian*)
			singlebodyZ=Table[ReplacePart[fullId,i->zOpList[[i]] ],{i,ncl}];
			singleHamiZ=Table[Subscript[Bz, i],{i,ncl}].Table[KroneckerProduct[##]&@@singlebodyZ[[i]],{i,ncl}];
			
			singlebodyX=Table[ReplacePart[fullId,i->xOpList[[i]] ],{i,ncl}];
			singleHamiX=Table[Subscript[Bx, i],{i,ncl}].Table[KroneckerProduct[##]&@@singlebodyX[[i]],{i,ncl}];
			
			singlebodyY=Table[ReplacePart[fullId,i->yOpList[[i]] ],{i,ncl}];
			singleHamiY=Table[Subscript[By, i],{i,ncl}].Table[KroneckerProduct[##]&@@singlebodyY[[i]],{i,ncl}];
			singleHami=singleHamiZ+singleHamiX+singleHamiY;
			
			
			(*two body  Hamiltonian*)
			xxRules={#[[1]]->xOpList[[ #[[1]] ]],#[[2]]->xOpList[[ #[[2]] ]]}&/@pairIdxlist;
			xxCoupling=Table[ReplacePart[fullId,xxRules[[i]]],{i,lpair}];
			xxHami=Table[Subscript[xx, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@xxCoupling[[i]],{i,lpair}];
			
			yyRules={#[[1]]->yOpList[[ #[[1]] ]],#[[2]]->yOpList[[ #[[2]] ]]}&/@pairIdxlist;
			yyCoupling=Table[ReplacePart[fullId,yyRules[[i]]],{i,lpair}];
			yyHami=Table[Subscript[yy, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@yyCoupling[[i]],{i,lpair}];
			
			zzRules={#[[1]]->zOpList[[ #[[1]] ]],#[[2]]->zOpList[[ #[[2]] ]]}&/@pairIdxlist;
			zzCoupling=Table[ReplacePart[fullId,zzRules[[i]]],{i,lpair}];
			zzHami=Table[Subscript[zz, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@zzCoupling[[i]],{i,lpair}];
			
			xyRules={#[[1]]->xOpList[[ #[[1]] ]],#[[2]]->yOpList[[ #[[2]] ]]}&/@pairIdxlist;
			xyCoupling=Table[ReplacePart[fullId,xyRules[[i]]],{i,lpair}];
			xyHami=Table[Subscript[xy, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@xyCoupling[[i]],{i,lpair}];
			
			
			yxRules={#[[1]]->yOpList[[ #[[1]] ]],#[[2]]->xOpList[[ #[[2]] ]]}&/@pairIdxlist;
			yxCoupling=Table[ReplacePart[fullId,yxRules[[i]]],{i,lpair}];
			yxHami=Table[Subscript[yx, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@yxCoupling[[i]],{i,lpair}];
			
			yzRules={#[[1]]->yOpList[[ #[[1]] ]],#[[2]]->zOpList[[ #[[2]] ]]}&/@pairIdxlist;
			yzCoupling=Table[ReplacePart[fullId,yzRules[[i]]],{i,lpair}];
			yzHami=Table[Subscript[yz, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@yzCoupling[[i]],{i,lpair}];
			
			zyRules={#[[1]]->zOpList[[ #[[1]] ]],#[[2]]->yOpList[[ #[[2]] ]]}&/@pairIdxlist;
			zyCoupling=Table[ReplacePart[fullId,zyRules[[i]]],{i,lpair}];
			zyHami=Table[Subscript[zy, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@zyCoupling[[i]],{i,lpair}];
			
			zxRules={#[[1]]->zOpList[[ #[[1]] ]],#[[2]]->xOpList[[ #[[2]] ]]}&/@pairIdxlist;
			zxCoupling=Table[ReplacePart[fullId,zxRules[[i]]],{i,lpair}];
			zxHami=Table[Subscript[zx, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@zxCoupling[[i]],{i,lpair}];
			
			xzRules={#[[1]]->xOpList[[ #[[1]] ]],#[[2]]->zOpList[[ #[[2]] ]]}&/@pairIdxlist;
			xzCoupling=Table[ReplacePart[fullId,xzRules[[i]]],{i,lpair}];
			xzHami=Table[Subscript[xz, pairIdxlist[[i]]],{i,lpair}].Table[KroneckerProduct[##]&@@xzCoupling[[i]],{i,lpair}];
			
			res=singleHami+xxHami+yyHami+zzHami+xyHami+yzHami+zxHami+yxHami+zyHami+xzHami;
		];
		Return[res];
	]/;Length[Jvallist]==ncl

NumericClusterHami[hflist_,dipList_,hamiform_]:=
	Module[{res,hfXlist,hfYlist,hfZlist,dipXXList,dipYYList,dipZZList,dipXYList,dipYZList,dipZXList,dipXZList,dipYXList,dipZYList,
			rulesSglbody,rulesTwobody},
		
		{hfXlist,hfYlist,hfZlist}=hflist;
		{dipXXList,dipYYList,dipZZList,
		dipXYList,dipYZList,dipZXList,
		dipYXList,dipZYList,dipXZList}=dipList;
		
		rulesSglbody={Subscript[Bx, i_]:>hfXlist[[i]],Subscript[By, i_]:>hfYlist[[i]],Subscript[Bz, i_]:>hfZlist[[i]]};
		rulesTwobody={
			Subscript[xx, {i_,j_}]:>dipXXList[[i,j]],
			Subscript[yy, {i_,j_}]:>dipYYList[[i,j]],
			Subscript[zz, {i_,j_}]:>dipZZList[[i,j]],
			Subscript[xy, {i_,j_}]:>dipXYList[[i,j]],
			Subscript[yz, {i_,j_}]:>dipYZList[[i,j]],
			Subscript[zx, {i_,j_}]:>dipZXList[[i,j]],
			Subscript[yx, {i_,j_}]:>dipYXList[[i,j]],
			Subscript[zy, {i_,j_}]:>dipZYList[[i,j]],
			Subscript[xz, {i_,j_}]:>dipXZList[[i,j]]
		};
		res=N[hamiform/.rulesSglbody/.rulesTwobody];
		Return[res];
	]

Options[ClusterHami]={MediatedTensor->False};
ClusterHami[coord1_,state_,b_,theta_,phi_,opt___]:=
	Module[{res,len,(*element,*)spinvalList,hamiform,coord,hfx,hfy,hfz,zeeman,bxList,byList,bzList,
		dipData,dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX,gammalist,
		mediatedintData,medXX,medYY,medZZ,medXY,medYZ,medZX,medYX,medZY,medXZ,
		isMed},
		
		isMed=MediatedTensor/.{opt}/.Options[ClusterHami];
		
		len=coord1//Length;

		{coord,spinvalList,gammalist}=CoordConvert[coord1];
		gammalist=gammalist/$GammaC13;
		hamiform=HamiNSpinXYZ[len,spinvalList];
		
		zeeman=Table[-b*$GammaC13*gammalist[[i]]/1000.0,{i,1,len}];
		{hfx,hfy,hfz}=Transpose[Table[gammalist[[i]]*HyperfineField[coord[[i]],state,b,theta,phi],{i,1,len}]];
		{bxList,byList,bzList}={hfx,hfy,hfz+zeeman};

		dipData=Table[
			If[i==j,
				{0.,0.,0., 0.,0.,0.}
				,
				gammalist[[i]]*gammalist[[j]]*DipolarInteraction[coord[[i]],coord[[j]],b,theta,phi]
			]
			,{i,1,len},{j,1,len}
		];
		
		If[isMed,
			mediatedintData=Table[
				If[i==j,
					{0.,0.,0., 0.,0.,0., 0.,0.,0.}
					,
					gammalist[[i]]*gammalist[[j]]*Flatten[MediatedInteraction[coord[[i]],coord[[j]],state,b,theta,phi]]
				]
				,{i,1,len},{j,1,len}
			];
			,
			mediatedintData=Table[{0.,0.,0., 0.,0.,0., 0.,0.,0.},{i,1,len},{j,1,len}];
		];
		
		{dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX}=Table[dipData[[All,All,i]],{i,1,6}];
		
		{medXX,medYY,medZZ,medXY,medYZ,medZX,medYX,medZY,medXZ}=Table[mediatedintData[[All,All,i]],{i,1,9}];
		
		res=NumericClusterHami[{bxList,byList,bzList},
			{dipXX,dipYY,dipZZ, dipXY,dipYZ,dipZX, dipXY,dipYZ,dipZX}
			+{medXX,medYY,medZZ, medXY,medYZ,medZX, medYX,medZY,medXZ}
			,hamiform];
		Return[res];
	](*/;Length[coord1]<=6*)

NoiseMat[coord1_,state_,b_,theta_,phi_]:=
	Module[{res,coord,(*element,*)spinvalList,gammalist,hamiform,len,hfx,hfy,hfz,dipData,dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX},
		len=coord1//Length;

		{coord,spinvalList,gammalist}=CoordConvert[coord1];
		gammalist=gammalist/$GammaC13;
		hamiform=HamiNSpinXYZ[len,spinvalList];

		{hfx,hfy,hfz}=Transpose[Table[gammalist[[i]]*HyperfineField[coord[[i]],state,b,theta,phi],{i,1,len}]];
		
		dipData=Table[{0.,0.,0.,0.,0.,0.},{i,1,len},{j,1,len}];
		{dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX}=Table[dipData[[All,All,i]],{i,1,6}];

		res=NumericClusterHami[{hfx,hfy,hfz},{dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX,dipXY,dipYZ,dipZX},hamiform];
		Return[res];
	]
	
BathMat[coord1_,b_,theta_,phi_]:=
	Module[{res,coord,(*element,*)spinvalList,gammalist,hamiform,len,hfx,hfy,hfz,dipData,dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX,
		bxList,byList,bzList,zeeman},
		len=coord1//Length;

		{coord,spinvalList,gammalist}=CoordConvert[coord1];
		gammalist=gammalist/$GammaC13;
		hamiform=HamiNSpinXYZ[len,spinvalList];

		zeeman=Table[-b*$GammaC13*gammalist[[i]]/1000.0,{i,1,len}];
		{hfx,hfy,hfz}=Transpose[Table[{0.,0.,0.},{i,1,len}]];
		{bxList,byList,bzList}={hfx,hfy,hfz+zeeman};
		
		dipData=Table[
			If[i==j,
				{0.,0.,0.,0.,0.,0.}
				,
				gammalist[[i]]*gammalist[[j]]*DipolarInteraction[coord[[i]],coord[[j]],b,theta,phi]
			]
			,{i,1,len},{j,1,len}
		];
		
		{dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX}=Table[dipData[[All,All,i]],{i,1,6}];
		
		res=NumericClusterHami[{bxList,byList,bzList},{dipXX,dipYY,dipZZ,dipXY,dipYZ,dipZX,dipXY,dipYZ,dipZX},hamiform];
		Return[res];
	]	

End[];

(*=================================================================*)
(*Cluster Evolution*)
Begin["`ClusterEvolution`"];

TimeIntervalSeq[seq_]:=
	Module[{res,intervals,doubleseq,seq1,seq2},
		intervals=Differences[seq];
		doubleseq=(-intervals)~Join~Reverse[intervals];
		{seq1,seq2}=Partition[doubleseq,2]//Transpose;
		res={seq1,seq2};
		Return[res];
	]
	
matOp[mat_,para_]:=
	Module[{res,expval,wmat},
		{expval,wmat}=para;
		res=wmat.(expval*mat);
		Return[res];
	]

coherenceATt[n_,wlist_,expValsLeft_,expValsRight_,idmat_]:=
	Module[{res,expVals,foldpara,resmat},
		expVals=Riffle[expValsLeft,expValsRight];
		foldpara=Transpose[{expVals,wlist}];
		(*resmat=Fold[matOp,idmat,foldpara];*)
		resmat=Fold[#2[[2]].(#2[[1]]*#1)&,idmat,foldpara];
		
		res=Tr[resmat];
		Return[res];
	]

ClusterCoherence[npulse_,seq_,hami1_,hami2_,{tmin_,tmax_,dt_},func_:Abs]:=
	Module[{res,dim,e1,v1,e2,v2,Wmat,WmatList,intervalLeft,intervalRight,timelist,expValLeft,expValRight,ntime,
			resCoh,wlist,wlistTrans,idmat},
		dim=hami1//Length;
		idmat=N[IdentityMatrix[dim]/dim];
		{intervalLeft,intervalRight}=TimeIntervalSeq[seq];
		{e1,v1}=Eigensystem[hami1]//Chop;
		{e2,v2}=Eigensystem[hami2]//Chop;

		Wmat=Conjugate[v1].Transpose[v2];
		wlist=Table[Wmat,{npulse+1}];
		wlistTrans=Table[Wmat//ConjugateTranspose,{npulse+1}];
		WmatList=Riffle[wlistTrans,wlist];

		timelist=Table[t,{t,tmin,tmax,dt}];
		expValLeft=Table[Exp[-I*t*#*e2]&/@intervalLeft,{t,tmin,tmax,dt}];
		expValRight=Table[Exp[-I*t*#*e1]&/@intervalRight,{t,tmin,tmax,dt}];

		ntime=timelist//Length;

		resCoh=Table[coherenceATt[npulse,WmatList,expValLeft[[i]],expValRight[[i]],idmat],{i,1,ntime}];
		res=Transpose[{timelist,resCoh//func}];
		Return[res];
	]

(*Mathematica realization*)
Options[CoherentClusterEvolution]={MediatedTensor->False,PostFunction->Abs};
CoherentClusterEvolution[nspin_,coord1_,{trans1_,trans2_},{b_,theta_,phi_},
		seq_,{tmin_,tmax_,dt_},opt___]:=
	Module[{res,hami1,hami2,npulse,isMed,func},
		func=PostFunction/.{opt}/.Options[CoherentClusterEvolution];
		isMed=MediatedTensor/.{opt}/.Options[CoherentClusterEvolution];
		
		hami1=ClusterHami[coord1,trans1,b,theta,phi,MediatedTensor->isMed];
		hami2=ClusterHami[coord1,trans2,b,theta,phi,MediatedTensor->isMed];
		npulse=Length[seq]-2;
		
		res=ClusterCoherence[npulse,seq,hami1,hami2,{tmin,tmax,dt},func];
		Return[res];
	]

(*MathLink realization*)
CoherentClusterEvolution[nspin_,coord1_,{trans1_,trans2_},{b_,theta_,phi_},
	pulsename_,npulse_,{tmin_,tmax_,dt_},isBlk_]:=
	Module[{res,coord,lk,pulsetype,timelist,spinlist,gammalist,data},
		pulsetype=Which[
			pulsename=="UDD",
			1,
			
			pulsename=="CPMG",
			2,
			
			True,(*UDD by default*)
			1
		];
		timelist=Range[tmin,tmax,dt];
		
		{coord,spinlist,gammalist}=CoordConvert[coord1];
		spinlist=spinlist*2+1;
		
		lk=Install[$MRMathLinkPath<>"/MLClusterEvolution/MLClusterEvolution"];
		data=ClusterEvolution[coord//Flatten, spinlist, gammalist, {b,theta,phi}, timelist, npulse, pulsetype, {trans1, trans2}, isBlk];
		Uninstall[lk];
		
		res=Transpose[{timelist,data}];
		Return[res];
	]/; ( isBlk==0 || isBlk == 1 )

FilterFunc[w_,t_,seq_]:=Abs[Sum[(-1)^k (Exp[I w t*seq[[k+2]]]-Exp[I w t*seq[[k+1]]]),{k,0,(seq//Length)-2}]]^2;

FilterFunctionArray[wt_,seq_]:=
	Module[{res,lk},
		lk=Install[$MRMathLinkPath<>"FilterFunction/FilterFunction"];
		res=FilterFunction[#, 1., seq]&/@wt;
		Uninstall[lk];
		Return[res];
	]

End[];


Begin["`GaussianApproximation`"];
(*
SpectWeight[eigVec_,i_,j_,noisemat_]:=
	Module[{res,vL,vR,matele,dim},
		dim=eigVec//Length;
		vL=eigVec[[i]];
		vR=eigVec[[j]];
		matele=Conjugate[vL].noisemat.vR;
		res=4.0*Abs[matele]^2/dim;
		Return[res];
	]

(*NoiseWeightFreq[nspin_,coord_,{trans1_,trans2_},{b_,theta_,phi_},seq_,{tmin_,tmax_,dt_},zerostate_]:=*)
NoiseWeightFreq[coord1_,{trans1_,trans2_},{b_,theta_,phi_}]:=
	Module[{res,nspin,dim,HBhami,noisematrix1,noisematrix2,hb,bmat,eigenVal,eigenVec,freqweightList},
		nspin = coord1//Length;
		dim=2^nspin;
		
		HBhami=BathMat[coord1,b,theta,phi];
		noisematrix1=NoiseMat[coord1,trans1,b,theta,phi];
		noisematrix2=NoiseMat[coord1,trans2,b,theta,phi];
		
		hb=HBhami+0.5*(noisematrix1+noisematrix2);
		bmat=0.5*(noisematrix1-noisematrix2);
		
		{eigenVal,eigenVec}=Eigensystem[hb]//Chop;
		eigenVec=eigenVec[[Ordering[eigenVal]]];
		eigenVal=Sort[eigenVal];
		freqweightList=Flatten[Table[Flatten[{i,j,eigenVal[[j]]-eigenVal[[i]],SpectWeight[eigenVec,i,j,bmat]}],{i,1,dim},{j,i+1,dim}],1];
		res=Sort[freqweightList,Abs[#1[[3]]]<Abs[#2[[3]]]&];
		Return[res];
	](*/;nspin>=4*)
*)
NoiseAnalysis[coord1_,{trans1_,trans2_},{b_,theta_,phi_}]:=
	Module[{res,nspin,dim,HBhami,noisematrix1,noisematrix2,hb,bmat,eigenVal,eigenVec},
		nspin = coord1//Length;
	
		HBhami=BathMat[coord1,b,theta,phi];
		noisematrix1=NoiseMat[coord1,trans1,b,theta,phi];
		noisematrix2=NoiseMat[coord1,trans2,b,theta,phi];
		
		dim=HBhami//Length;
		
		hb=HBhami+0.5*(noisematrix1+noisematrix2);
		bmat=noisematrix1-noisematrix2;
		
		{eigenVal,eigenVec}=Eigensystem[hb]//Chop;
		eigenVec=eigenVec[[Ordering[eigenVal]]];
		eigenVal=Sort[eigenVal];
		
		res=Table[NoiseFromState[i,eigenVal,eigenVec,bmat],{i,1,dim}];
	
		Return[res];
	]

NoiseFromState[init_,val_,vec_,noisemat_] := 
	Module[{res,dim,vl,vr,ele2,freq,weight}, 
		dim = val//Length;
		
		res=Table[
			vl=vec[[init]];
			vr=vec[[i]];
			ele2=Abs[Conjugate[vl].noisemat.vr]^2;
			freq=Chop[val[[i]]-val[[init]]];
			weight=If[Abs[freq]>0.0, 0.5 * ele2/freq^2, Infinity];
			{init, i, freq, ele2, weight}
			,{i,1,dim}];
		
		Return[res];
	]
evolutionFromState[specinfo_,seq_,timelist_] := 
	Module[{res,dim,freqWeightPair,len,t,exSum,lk}, 
		dim=specinfo//Length;
		freqWeightPair=Select[specinfo[[All,{3,5}]], NumberQ[#[[2]]] &];
		
		len=timelist//Length;

		lk=Install[$MRMathLinkPath<>"FilterFunction/FilterFunction"];
		res=Table[
				t=timelist[[i]];
				exSum=Plus@@( (#[[2]] * FilterFunction[ #[[1]], t, seq ])&/@freqWeightPair );
				Exp[-exSum]
			,{i,1,len}];
		Uninstall[lk];
		Return[res];
	]
	
CoherentClusterEvolutionGaussian[nspin_,coord1_,{trans1_,trans2_},{b_,theta_,phi_},
		seq_,{tmin_,tmax_,dt_}] := 
	Module[{res, specInfo, dim, resInit,timelist},
		specInfo=NoiseAnalysis[coord1,{trans1,trans2},{b,theta,phi}];
		dim=specInfo//Length;
		timelist=Range[tmin,tmax,dt];
		resInit=Table[
			evolutionFromState[specInfo[[i]],seq,timelist],
			{i,1,dim}];
			
		res=Transpose[{timelist,Mean[resInit]}];
		
		Return[res];		
	]

CoherentClusterEvolutionGaussian1[nspin_,coord1_,{trans1_,trans2_},{b_,theta_,phi_},
		seq_,{tmin_,tmax_,dt_}] := 
	Module[{res,dim,HBhami,noisematrix1,noisematrix2,hb,bmat,eigenVal,eigenVec,freqweightPair,resExp,timelist,ele2,lk},
			
		HBhami=BathMat[coord1,b,theta,phi];
		noisematrix1=NoiseMat[coord1,trans1,b,theta,phi];
		noisematrix2=NoiseMat[coord1,trans2,b,theta,phi];
		
		dim=HBhami//Length;
		
		hb=HBhami+0.5*(noisematrix1+noisematrix2);
		bmat=noisematrix1-noisematrix2;
		
		{eigenVal,eigenVec}=Eigensystem[hb]//Chop;
		eigenVec=eigenVec[[Ordering[eigenVal]]];
		eigenVal=Sort[eigenVal];
		
		freqweightPair=Flatten[Table[
			ele2=Abs[Conjugate[ eigenVec[[i]] ].bmat.eigenVec[[j]] ]^2;
			{eigenVal[[j]]-eigenVal[[i]], ele2}
			,{i,1,dim},{j,i+1,dim}],1];
		freqweightPair=Select[freqweightPair,Abs[ #[[1]] ] >1.0*10^(-10)&];
						
		timelist=Range[tmin,tmax,dt];
		lk=Install[$MRMathLinkPath<>"FilterFunction/FilterFunction"];
		resExp=Table[
			Plus@@( ( #[[2]] / #[[1]]^2 * FilterFunction[ #[[1]], timelist[[i]], seq ]) &/@freqweightPair)
			,{i,1,timelist//Length}];
		Uninstall[lk];
			
		res=Transpose[{timelist,Exp[- resExp/dim ]}];
		
		Return[res];		
	]

(*
GaussianCoherentClusterEvolution[nspin_,{{i0_,j0_,k0_},r0_,d0_},{trans1_,trans2_},{b_,theta_,phi_},
		seq_,{tmin_,tmax_,dt_},zerostate_,Sel_,gammaBath_:$GammaC13]:=
	Module[{res,freqWeightPair,SelfreqWeightPair,len,freq,weight,expo},

		freqWeightPair=NoiseWeightFreq[nspin,{{i0,j0,k0},r0,d0},{trans1,trans2},{b,theta,phi},seq,{tmin,tmax,dt},zerostate];
		SelfreqWeightPair=freqWeightPair[[Sel]];
		len=SelfreqWeightPair//Length;
		
		expo=Table[
			Sum[
				{freq,weight}=SelfreqWeightPair[[i,{2,3}]];
				weight*FilterFunction[freq//Re,t,seq]
				,{i,1,len}
				]
			,{t,tmin,tmax,dt}
			];
		
		res=Exp[-expo];
		Return[res];
	]/;nspin==2||nspin==3

GaussianCoherentClusterEvolution[nspin_,coord_,{trans1_,trans2_},{b_,theta_,phi_},
		seq_,{tmin_,tmax_,dt_},zerostate_,Sel_,gammaBath_:$GammaC13]:=
	Module[{res,freqWeightPair,SelfreqWeightPair,len,freq,weight,expo},

		freqWeightPair=NoiseWeightFreq[nspin,coord,{trans1,trans2},{b,theta,phi},seq,{tmin,tmax,dt},zerostate];
		SelfreqWeightPair=freqWeightPair[[Sel]];
		len=SelfreqWeightPair//Length;
		
		expo=Table[
			Sum[
				{freq,weight}=SelfreqWeightPair[[i,{2,3}]];
				weight*FilterFunction[freq//Re,t,seq]
				,{i,1,len}
				]
			,{t,tmin,tmax,dt}
			];
		
		res=Exp[-expo];
		Return[res];
	]/;nspin>=4
*)
End[];

(*=================================================================*)
(*Load CCE Results*)
Begin["`LoadCCE`"];

ReadCCEData[filename_]:=
	Module[{res,st,ntime,cceOrder,clstNum,singNum,clstIndex,nspin,spinIdx,redval,irredval,clsInfo},

		st=OpenRead[filename,BinaryFormat->True];
		{ntime,cceOrder}=BinaryRead[st,{"Integer32","Integer32"}];
		clstNum=BinaryRead[st,Table["Integer32",{cceOrder}]];
		singNum=BinaryRead[st,"Integer32"];

		clsInfo=Table[
			{clstIndex,nspin}=BinaryRead[st,{"Integer32","Integer32"}];
			spinIdx=BinaryRead[st,Table["Integer32",{nspin}]];
			redval=BinaryRead[st,Table["Real64",{ntime+1}]];
			irredval=BinaryRead[st,Table["Real64",{ntime+1}]];
			{clstIndex,nspin,spinIdx,redval,irredval}
			,{i,1,singNum}
		];
		
		Close[st];
		
		res={ntime,cceOrder,clstNum,singNum,clsInfo};
		Return[res];
	]

PlotCCECluster[timelist_,ccedata_,i_,isRed_,idxShift_,opt___]:=
	Module[{res,data,reducible,pts,label},
		reducible=If[isRed=="R",4,5];
		data=ccedata[[-1,i,reducible]];
		pts=Transpose[{timelist,data}];
		label=ccedata[[-1,i,{1,2,3}]]+{0,0,idxShift};
		res=ListPlot[pts,PlotLabel->ToString[label],opt];
		Return[res];
	];

GlobalCCEInfo[ccedata_]:=
	Module[{res,cceorder,names,data},
		cceorder=ccedata[[2]];
		names={"NTime","CCE Order"}~Join~Table["CCE-"<>ToString[i],{i,cceorder}]~Join~{"Singular Num"};
		data=ccedata[[1;;4]]//Flatten;
		res=Transpose[{names,data}];
		Return[res];
	]

End[];

(*=================================================================*)
EndPackage[];
